{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "03995d6e-9a58-4c26-9b9a-a789cb90dd90",
   "metadata": {},
   "source": [
    "# Using QSVT for fixed-point amplitude amplification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdf4c557-fbf3-43ad-a840-4d25a038931e",
   "metadata": {},
   "source": [
    "This demo will show how to use the QSVT framework for search problems. Specifically, we will implement fixed-point amplitude amplification (FPAA). In FPAA, we do not know in advance the concentration of solutions to the search problem, want to sample a solution with high probability. In contrast, for the original grover search algorithm, too much iterations might 'overshoot'.\n",
    "\n",
    "The demo is based on the paper [Grand unification of quantum algorithms](#grand).\n",
    "\n",
    "Given $|s\\rangle$ the initial state and $|t\\rangle$ the 'good' states, we get an effective block encoding of a 1-dimensional matrix $A=|t\\rangle\\langle s|$.\n",
    "\n",
    "Given that  $a = \\langle s|t\\rangle\\gt0$, we want to amplify $a$. The signal operator $U$ here will be $I$ (and also $\\dagger{U}$). Now we implement 2 projector-rotations - one in '$|s\\rangle$' space and one in '$|t\\rangle$' space, each one around the given state, giving phase to the specific state."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "170f965a-4bf2-4a4a-8734-03e4d1c6e0b0",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Defining the QSVT circuit for the problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e456f39b-39a5-4232-9bbf-fbc6224dfad7",
   "metadata": {},
   "source": [
    "We start with the general qsvt framework definition. It accepts a unitary which block-encode a matrix together with projector-controlled-phase functions which rotate the state around each of the subspaces where the matrix is encoded.\n",
    "\n",
    "It applies the `qsvt_step` multiple times, iterating over the rotation angles provided which encode the polynomial transformation.\n",
    "Notice - The last step is quite tricky and depend on the specific transformation we wish to perform. Here the code is suitable for the FP Amplitude Amplification case. Also - we wrap the auxilliary qubit with $H$ gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "352c9910-50ed-4560-b371-a94992ca76b0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:33.089927Z",
     "iopub.status.busy": "2024-05-07T14:31:33.087719Z",
     "iopub.status.idle": "2024-05-07T14:31:35.868598Z",
     "shell.execute_reply": "2024-05-07T14:31:35.867923Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import QArray, QBit, QNum, qfunc\n",
    "from classiq.qmod import (\n",
    "    IDENTITY,\n",
    "    RZ,\n",
    "    CArray,\n",
    "    CInt,\n",
    "    CReal,\n",
    "    H,\n",
    "    QCallable,\n",
    "    X,\n",
    "    allocate,\n",
    "    apply_to_all,\n",
    "    control,\n",
    "    hadamard_transform,\n",
    "    if_,\n",
    "    invert,\n",
    "    repeat,\n",
    ")\n",
    "from classiq.qmod.symbolic import logical_and\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def my_qsvt_step(\n",
    "    phase_seq: CArray[CReal],\n",
    "    index: CInt,\n",
    "    proj_ctrl_phase_1: QCallable[CReal, QArray[QBit], QBit],\n",
    "    proj_ctrl_phase_2: QCallable[CReal, QArray[QBit], QBit],\n",
    "    u: QCallable[QArray[QBit]],\n",
    "    qbv_reg: QArray[QBit],\n",
    "    qbv_aux: QBit,\n",
    "):\n",
    "    proj_ctrl_phase_1(phase_seq[2 * index], qbv_reg, qbv_aux)\n",
    "    u(qbv_reg)\n",
    "    proj_ctrl_phase_2(phase_seq[2 * index + 1], qbv_reg, qbv_aux)\n",
    "\n",
    "    if_(\n",
    "        condition=2 * index + 2 == phase_seq.len,\n",
    "        then=lambda: IDENTITY(qbv_reg),\n",
    "        else_=lambda: invert(lambda: u(qbv_reg)),\n",
    "    )\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def my_qsvt(\n",
    "    phase_seq: CArray[CReal],\n",
    "    proj_ctrl_phase_1: QCallable[CReal, QArray[QBit], QBit],\n",
    "    proj_ctrl_phase_2: QCallable[CReal, QArray[QBit], QBit],\n",
    "    u: QCallable[QArray[QBit]],\n",
    "    qbv_reg: QArray[QBit],\n",
    "    qbv_aux: QBit,\n",
    ") -> None:\n",
    "    H(qbv_aux)\n",
    "\n",
    "    repeat(\n",
    "        count=phase_seq.len / 2,\n",
    "        iteration=lambda index: my_qsvt_step(\n",
    "            phase_seq, index, proj_ctrl_phase_1, proj_ctrl_phase_2, u, qbv_reg, qbv_aux\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    H(qbv_aux)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba2ae29c-e041-45a3-abf2-be7da1814650",
   "metadata": {},
   "source": [
    "## Initial state z-rotation - rotation around $|s\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2066e586-3b7e-4a4c-b1eb-42499026a509",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:35.871793Z",
     "iopub.status.busy": "2024-05-07T14:31:35.871221Z",
     "iopub.status.idle": "2024-05-07T14:31:35.875100Z",
     "shell.execute_reply": "2024-05-07T14:31:35.874586Z"
    }
   },
   "outputs": [],
   "source": [
    "def initial_state_rot(phase: CReal, state: QArray[QBit], aux: QBit):\n",
    "    hadamard_transform(state)\n",
    "    apply_to_all(X, state)\n",
    "    control(ctrl=state, operand=lambda: X(aux))\n",
    "    RZ(phase, aux)\n",
    "    control(ctrl=state, operand=lambda: X(aux))\n",
    "    apply_to_all(X, state)\n",
    "    hadamard_transform(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f93e8ef7-cc2f-423b-9bbe-8365c7d0b3ba",
   "metadata": {},
   "source": [
    "## Good states z-rotation - rotation around $|t\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8fe5fdfa-b63a-4a5b-a6db-d947da2e6f63",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:35.877471Z",
     "iopub.status.busy": "2024-05-07T14:31:35.877060Z",
     "iopub.status.idle": "2024-05-07T14:31:35.880346Z",
     "shell.execute_reply": "2024-05-07T14:31:35.879911Z"
    }
   },
   "outputs": [],
   "source": [
    "def target_state_rot(\n",
    "    phase: CReal,\n",
    "    arith_oracle: QCallable[QArray[QBit], QBit],\n",
    "    state: QArray[QBit],\n",
    "    aux: QBit,\n",
    "):\n",
    "    arith_oracle(state, aux)\n",
    "    RZ(phase, aux)\n",
    "    arith_oracle(state, aux)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40488ffb-ebd4-49b4-aa24-7ae212150202",
   "metadata": {},
   "source": [
    "## Defining the arithmetic Oracle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dd55e04b-443d-409c-bd3a-8e57aeedef95",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:35.883461Z",
     "iopub.status.busy": "2024-05-07T14:31:35.883164Z",
     "iopub.status.idle": "2024-05-07T14:31:35.892579Z",
     "shell.execute_reply": "2024-05-07T14:31:35.891639Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "EQUATION = \"(a + b) == 3 and (c - a) == 2\"\n",
    "\n",
    "REG_SIZES = {\"a\": 2, \"b\": 1, \"c\": 3}\n",
    "REG_SIZE = sum(REG_SIZES.values())\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def arith_equation(a: QNum, b: QNum, c: QNum, res: QBit):\n",
    "    res ^= logical_and((a + b) == 3, (c - a) == 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d727b48b-672d-4734-8b98-348943ce4fb6",
   "metadata": {},
   "source": [
    "## Wrapping everything together for the FPAA case"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e79829-8e88-448e-824b-f1bc72470d81",
   "metadata": {},
   "source": [
    "In the FPAA case, the provided unitary is just the Identity matrix! In addition, we provide both projector-controlled phase functions - the initial and target state rotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "863fda93-7371-4265-a02e-8da9939e659e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:35.897345Z",
     "iopub.status.busy": "2024-05-07T14:31:35.896211Z",
     "iopub.status.idle": "2024-05-07T14:31:35.904983Z",
     "shell.execute_reply": "2024-05-07T14:31:35.904370Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq.qmod import Output\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def qsvt_fpaa(\n",
    "    phase_seq: CArray[CReal],\n",
    "    arith_oracle: QCallable[QArray[QBit], QBit],\n",
    "    qbv_reg: QArray[QBit],\n",
    "    qbv_aux: Output[QBit],\n",
    ") -> None:\n",
    "    allocate(1, qbv_aux)\n",
    "\n",
    "    my_qsvt(\n",
    "        phase_seq,\n",
    "        lambda phase, state, aux: initial_state_rot(phase, state, aux),\n",
    "        lambda phase, state, aux: target_state_rot(phase, arith_oracle, state, aux),\n",
    "        lambda state: IDENTITY(target=state),\n",
    "        qbv_reg,\n",
    "        qbv_aux,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "640e9b8d-90eb-4caa-83d4-6f87be59a15f",
   "metadata": {},
   "source": [
    "## get the phase sequence for the sign function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62190792-6bca-434e-b5cf-e0a03cd4aac0",
   "metadata": {
    "tags": []
   },
   "source": [
    "Now we will use the package `pyqsp` in order to get the phases of the rotation sequence.\n",
    "Get directly the coef of the sign function, based on the erfc approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "73d83af6-8d2c-4e8a-b392-1ba337f74ca0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:35.909571Z",
     "iopub.status.busy": "2024-05-07T14:31:35.909192Z",
     "iopub.status.idle": "2024-05-07T14:31:35.930875Z",
     "shell.execute_reply": "2024-05-07T14:31:35.929879Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[pyqsp.poly.PolySign] degree=25, delta=10\n",
      "[PolyTaylorSeries] max [0.89999873] is at [1.]: normalizing\n",
      "[PolyTaylorSeries] average error = 0.0962788793562509 in the domain [-1, 1] using degree 25\n"
     ]
    }
   ],
   "source": [
    "import pyqsp\n",
    "\n",
    "DEGREE = 25\n",
    "X_BASIS = True\n",
    "\n",
    "pg = pyqsp.poly.PolySign()\n",
    "pcoefs, scale = pg.generate(\n",
    "    degree=DEGREE, delta=10, ensure_bounded=True, return_scale=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "267d3e27-135b-4167-b6b6-8a06a8a99aaa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:35.935069Z",
     "iopub.status.busy": "2024-05-07T14:31:35.933978Z",
     "iopub.status.idle": "2024-05-07T14:31:36.300678Z",
     "shell.execute_reply": "2024-05-07T14:31:36.299925Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAHACAYAAACyDu7BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABN+UlEQVR4nO3deXgUZb728bs7eyALELJBJEAQ1CAgynZAYMjIpgPCOIgzsgi4HGWOgKLoCG4jqIiio6IooOeIKBxZRhEEBBWNIDvIInvYAiaRhCSQpNP1/sHbfdImQKfpTqXD93NddaW7+qnqX1XK5s7j009ZDMMwBAAAAOCirGYXAAAAAPgDgjMAAADgBoIzAAAA4AaCMwAAAOAGgjMAAADgBoIzAAAA4AaCMwAAAOAGgjMAAADghkCzC6ju7Ha7jh8/roiICFksFrPLAQAAwO8YhqEzZ84oMTFRVqvv+oUJzpdw/PhxJSUlmV0GAAAALuHIkSNq2LChz/ZPcL6EiIgISed/EZGRkSZXAwAAgN/Ly8tTUlKSM7f5CsH5EhzDMyIjIwnOAAAA1Zivh9Xy5UAAAADADQRnAAAAwA0EZwAAAMANjHEGAAB+xzAM2Ww2lZaWml0KqkBAQIACAwNNnxqY4AwAAPxKcXGxTpw4ocLCQrNLQRUKDw9XQkKCgoODTauB4AwAAPyG3W7XwYMHFRAQoMTERAUHB5veCwnfMgxDxcXF+vXXX3Xw4EE1a9bMpzc5uRiCMwAA8BvFxcWy2+1KSkpSeHi42eWgioSFhSkoKEiHDx9WcXGxQkNDTamDLwcCAAC/Y1aPI8xTHX7n5lcAAAAA+AGCMwAAgB9Ys2aNLBaLLBaL+vfv77X9zpkzR9HR0ZXe7qmnntK9997rdvusrCzFxsbq6NGjLuuHDRvmPK5FixZVuo6qRHAGAACoAmUDYlBQkBo3bqzx48fr3LlzldrPnj17NGfOnAr3GxwcrJSUFD377LOy2WxePoL/k5mZqenTp+vJJ590e5uYmBgNGTJEkyZNclk/ffp0nThxwtsl+gTBGQAAoIr06tVLJ06c0IEDB/Tqq6/qnXfeKRckLyU2NrZcD7Fjv3v37tW4ceP09NNP6+WXX/Zi5a7ee+89derUSY0aNarUdsOHD9dHH32knJwc57qoqCjFx8d7u0SfIDgDAABUkZCQEMXHxyspKUn9+/dXWlqaVqxYIen8VHuTJ09W48aNFRYWplatWmnBggWV2m+jRo30wAMPKC0tTUuWLJEk/fbbbxoyZIjq1Kmj8PBw9e7dW3v37q1wP4cOHZLVatWGDRtc1r/22mtq1KiR7Ha7JGnevHm67bbbXNosW7ZMnTt3VnR0tOrVq6dbb71V+/fvd2lz3XXXKTExUQsXLnTruKobgjMAAPBrhmGooKDAlMUwDI/r3rFjh3744QfnDT0mT56sDz/8UDNmzNDPP/+sMWPG6G9/+5u++eabSu87LCxMxcXFks4P5diwYYOWLFmi9PR0GYahPn36qKSkpNx2ycnJSktL0+zZs13Wz549W8OGDZPValVOTo527typG2+80aVNQUGBxo4dqw0bNmjVqlWyWq26/fbbnWHboV27dvruu+8qfUzVAfM4AwAAv1ZYWKjatWub8t75+fmqVauW2+0///xz1a5dWzabTUVFRbJarfrXv/6loqIivfDCC1q5cqU6duwoSWrSpInWrl2rd955R127dnVr/4ZhaNWqVVq+fLlGjx6tvXv3asmSJfr+++/VqVMnSdJHH32kpKQkLVq0SHfccUe5fYwcOVL333+/pk2bppCQEG3atEnbt2/X4sWLJUkZGRkyDEOJiYku2w0cONDl+axZs1S/fn3t3LlTqampzvWJiYnavHmz2+esOqHHGQAAoIp0795dW7Zs0bp16zR06FANHz5cAwcO1L59+1RYWKg//vGPql27tnP58MMPyw13qIgjkIeGhqp3794aNGiQnn76ae3atUuBgYFq3769s229evXUvHlz7dq1q8J99e/fXwEBAc7hFHPmzFH37t2VnJwsSTp79qwklbsJyd69ezV48GA1adJEkZGRzvYZGRku7cLCwvz2dun0OAMAAL8WHh6u/Px80967MmrVqqWUlBRJ53tkW7Vqpffff9/ZI/vFF1+oQYMGLtuEhIRccr/du3fX22+/reDgYCUmJiow0POIFxwcrCFDhmj27NkaMGCA5s6dq+nTpztfj4mJkXR+7HT9+vWd62+77TY1atRIM2fOVGJioux2u1JTU51DRhxycnJctvMnBGcAAODXLBZLpYZLVBdWq1VPPPGExo4dq19++UUhISHKyMhwe1hGWWUDeVnXXHONbDab1q1b5xyqkZ2drT179ujaa6+94P5Gjhyp1NRUvfXWW7LZbBowYIDztaZNmyoyMlI7d+7U1Vdf7bLPmTNnqkuXLpKktWvXVrjvHTt2qFu3bpU+xuqAoRoAAAAmueOOOxQQEKB33nlHjzzyiMaMGaMPPvhA+/fv16ZNm/TGG2/ogw8+8Hj/zZo1U79+/TRq1CitXbtWW7du1d/+9jc1aNBA/fr1u+B211xzjTp06KDHHntMgwcPVlhYmPM1q9WqtLQ0l2Bcp04d1atXT++++6727dunr7/+WmPHji2338LCQm3cuFG33HKLx8dkJoIzAACASQIDA/XQQw/ppZde0oQJE/TUU09p8uTJuuaaa9SrVy998cUXaty48WW9x+zZs9W2bVvdeuut6tixowzD0NKlSxUUFHTR7UaMGKHi4mLdc8895V4bOXKk5s2b55wxw2q1at68edq4caNSU1M1ZsyYCueRXrx4sa666ipnr7S/sRiXM4/KFSAvL09RUVHKzc1VZGSk2eUAAHBFO3funA4ePKjGjRuX+3JaTbdmzRp1795dv/32m0e3yK6s5557TvPnz9e2bdvKvWYYhtq3b68xY8Zo8ODBbu+zQ4cO+vvf/6677rqr3GsWi0ULFy684O3EL/a7r6q8Ro8zAACAH2nYsGGlwmpl5efna8eOHfrXv/6l0aNHV9jGYrHo3XffrdRtvbOysjRgwIBytd9///2mTSdYWfQ4XwI9zgAAVB9Xco/z2bNndezYMUlS7dq1fXab6mHDhunjjz9W//79NXfuXAUEBPjkfRxOnTqlvLw8SVJCQsIFv+hZHXqcmVUDAADAD4SFhVU4c4a3zZkzR3PmzPH5+zjExsYqNja2yt7vcjBUAwAAAHADwRkAAABwA8EZAAAAcAPBGQAAAHADwRkAAABwA8EZAAAAcAPBGQAAwA+sWbNGFotFFovlgnfX88ScOXM8uhPhU089pXvvvdft9llZWYqNjdXRo0dd1g8bNsx5XIsWLap0HVWJ4AwAAFAFhg0b5pXAu2fPHpd5lssGz+DgYKWkpOjZZ5+t1F39KiszM1PTp0/Xk08+6fY2MTExGjJkiCZNmuSyfvr06Tpx4oS3S/QJgjMAAIAfiY2NLddD3KtXL504cUJ79+7VuHHj9PTTT+vll1/2WQ3vvfeeOnXqpEaNGlVqu+HDh+ujjz5STk6Oc11UVJTP7oLobQRnAACAKtatWzeNHj1aDz/8sOrUqaO4uDjNnDlTBQUFGj58uCIiIpSSkqIvv/zSrf2FhIQoPj5ejRo10gMPPKC0tDQtWbJEkvTbb79pyJAhqlOnjsLDw9W7d2/t3bu3wv0cOnRIVqtVGzZscFn/2muvqVGjRrLb7ZKkefPm6bbbbnNps2zZMnXu3FnR0dGqV6+ebr31Vu3fv9+lzXXXXafExEQtXLjQreOqbgjOAADArxmGVFBgzmIYntf9wQcfKCYmRuvXr9fo0aP1wAMP6I477lCnTp20adMm3XLLLbr77rtVWFhY6X2HhYWpuLhY0vmhHBs2bNCSJUuUnp4uwzDUp08flZSUlNsuOTlZaWlpmj17tsv62bNna9iwYbJarcrJydHOnTt14403urQpKCjQ2LFjtWHDBq1atUpWq1W33367M2w7tGvXTt99912lj6k6CDS7AAAAgMtRWCjVrm3Oe+fnS7VqebZtq1at9I9//EOSNGHCBE2ZMkUxMTEaNWqUJGnixIl6++23tW3bNnXo0MGtfRqGoVWrVmn58uUaPXq09u7dqyVLluj7779Xp06dJEkfffSRkpKStGjRIt1xxx3l9jFy5Ejdf//9mjZtmkJCQrRp0yZt375dixcvliRlZGTIMAwlJia6bDdw4ECX57NmzVL9+vW1c+dOpaamOtcnJiZq8+bNbp6l6oUeZwAAABNcf/31zscBAQGqV6+eWrZs6VwXFxcnSTp16tQl9/X555+rdu3aCg0NVe/evTVo0CA9/fTT2rVrlwIDA9W+fXtn23r16ql58+batWtXhfvq37+/AgICnMMp5syZo+7duys5OVmSdPbsWUlSaGioy3Z79+7V4MGD1aRJE0VGRjrbZ2RkuLQLCwvzqBe9OqDHGQAA+LXw8PM9v2a9t6eCgoJcnlssFpd1FotFksoNdahI9+7d9fbbbys4OFiJiYkKDPQ84gUHB2vIkCGaPXu2BgwYoLlz52r69OnO12NiYiSdHztdv3595/rbbrtNjRo10syZM5WYmCi73a7U1FTnkBGHnJwcl+38CcEZAAD4NYvF8+ESNUWtWrWUkpJSbv0111wjm82mdevWOYdqZGdna8+ePbr22msvuL+RI0cqNTVVb731lmw2mwYMGOB8rWnTpoqMjNTOnTt19dVXu+xz5syZ6tKliyRp7dq1Fe57x44d6tatm6eHaiqGagAAANRQzZo1U79+/TRq1CitXbtWW7du1d/+9jc1aNBA/fr1u+B211xzjTp06KDHHntMgwcPVlhYmPM1q9WqtLQ0l2Bcp04d1atXT++++6727dunr7/+WmPHji2338LCQm3cuFG33HKLdw+0ihCcAQAAarDZs2erbdu2uvXWW9WxY0cZhqGlS5eWGyryeyNGjFBxcbHuueeecq+NHDlS8+bNcw4jsVqtmjdvnjZu3KjU1FSNGTOmwnmkFy9erKuuusrZK+1vGKoBAABQBcre7W/NmjXlXj906FC5dYYb892V3W9F6tSpow8//PCCrw8bNkzDhg0rt/7YsWNq2bKlbrrppnKv9erVS4mJifrkk080ePBgSVJaWpp27tx50fqnT5+uiRMnXrTe6oweZwAAAD/SsGFDZ1j1hfz8fO3YsUP/+te/NHr06ArbWCwWvfvuu5W6rXdWVpYGDBhQrvb7779ftc2aT7CSLIY7f8pcwfLy8hQVFaXc3FxFRkaaXQ4AAFe0c+fO6eDBg2rcuHG56dBqurNnz+rYsWOSpNq1a/vsNtXDhg3Txx9/rP79+2vu3LkKCAjwyfs4nDp1Snl5eZKkhIQE1brANz0v9ruvqrzGUA0AAAA/EBYWVuHMGd42Z86cSw7/8KbY2FjFxsZW2ftdDr8aqvHtt9/qtttuU2JioiwWixYtWnTJbdasWaMbbrhBISEhSklJqdILAQAAADWHXwXngoICtWrVSm+++aZb7Q8ePKi+ffuqe/fu2rJlix5++GGNHDlSy5cv93GlAAAAqGn8aqhG79691bt3b7fbz5gxQ40bN9Yrr7wi6fychGvXrtWrr76qnj17+qpMAADgY3xF68pTHX7nftXjXFnp6elKS0tzWdezZ0+lp6dfcJuioiLl5eW5LAAAoHpwzD1cWFhociWoao7f+aXmn/Ylv+pxrqzMzEzFxcW5rIuLi1NeXp7Onj3rchcch8mTJ+uZZ56pqhIBAEAlBAQEKDo6WqdOnZIkhYeHy2KxmFwVfMkwDBUWFurUqVOKjo72+SwfF1Ojg7MnJkyY4HKLyLy8PCUlJZlYEQAAKMsxDZsjPOPKEB0d7bMp+NxVo4NzfHy8Tp486bLu5MmTioyMrLC3WZJCQkIUEhJSFeUBAAAPWCwWJSQkKDY2ViUlJWaXgyoQFBRkak+zQ40Ozh07dtTSpUtd1q1YsUIdO3Y0qSIAAOAtAQEB1SJM4crhV18OzM/P15YtW7RlyxZJ56eb27JlizIyMiSdH2YxZMgQZ/v7779fBw4c0Pjx47V792699dZb+vTTTzVmzBgzygcAAIAf86vgvGHDBrVp00Zt2rSRJI0dO1Zt2rTRxIkTJUknTpxwhmhJaty4sb744gutWLFCrVq10iuvvKL33nuPqegAAABQaRajOkyKV41V1b3PAQAA4Jmqymt+1eMMAAAAmIXgDAAAALiB4AwAAAC4oUZPRwcAuDKdOXNGBQUFKi4uVnFxsWrVqqV69eopODjY7NIA+DGCMwDAr2VmZmr58uX67rvvtHv3bu3Zs0dZWVkVtq1du7aSkpLUvHlzXX311c4lNTVVderUqeLKAfgbgjMAwO+UlpZq4cKFmjFjhlavXi273V6ujcViUUhIiAIDA1VYWCi73a78/Hzt2rVLu3btKte+YcOGuv7669WyZUvnz+bNm1ebXmqbzabMzEwdO3ZMR48e1bFjx3Ty5En99ttv5ZYzZ844e9uLi4udd9cLCgpScHCwy8/Q0FBFRUW5LNHR0S6PK1oiIyNNuflISUmJCgsLVVBQ4PKz7OOioiLZ7XaXpbS01Pk4ICCg3Hko+zM4ONh5J+ELLUFBQbJYLF4/PsMwVFpaKpvNppKSEpflQutsNpusVquCgoIUGBiooKAgl8e//xkQEKDAwEAFBATIanV/1G7Z2hw/y57b0tJSl8eO8y2d/+/x98uF1l+qreNx2Z9nzpzx2u/gYpiO7hKYjg4Aqg/DMDR//nw9+eST2rdvn3P9jTfeqLS0NLVq1UrNmzdXs2bNVLt2befrdrtdubm5ysrK0sGDB/XLL79oz549+uWXX7R7926XewCUFRQUpGuuucYZph2BOjEx0auhqaioSMePH9fRo0d19OhRHTlyxPnYEZIzMzMr/APBTBEREeUCdVRUlEJCQpzh0xFEg4ODncHr92Hr3Llz5cLvhR5Xp1tsXyxYS7poeL9YIK5qjhBdNlAbhuESkB21V3e+zmsE50sgOANA9XDgwAGNGDFCa9askSTVq1dPDzzwgO655x41btz4svZ9+vRp7dixQ9u2bdP27dudPy/UixUeHq7ExEQ1aNBADRo0UGJioiIjIxUWFqbw8HCFh4crKChIJSUlLj2/eXl5ys7Odi5ZWVnOnmN3BAYGKiEhQQ0bNlSDBg0UHx+vOnXqlFsc4bVsT6okZz1lf549e1a5ubkuy+nTp8s9L7sUFhZe1vn2BqvVqvDwcNWqVcvlZ3h4uEJCQpy9qVar1eWx1WpVaWlphefC8XsqKipSUVGRy+OioiLZbDbTjvf3vchll4CAANntdpcA7njs+FlaWurzGh3n2nG+AwICnH9gGobhslS07kJLZRCcTUZwBgDzffzxxxo1apQKCgoUHh6u8ePHa9y4cS69yt5mGIYOHz5cLkzv2bPHJz1vISEhSkpKUsOGDcstjoAeGxtryvCI3ysuLq4wUJ8+fVp5eXnO0Fk2iBYXF8tisZQLVlarVaGhoc7QWzYAl338+5AcHBzsk6ESF1NaWlouTF9osVgsLmG9bIi3WCzlwm9Fgbjs0IrLPVZHj7cjRJftSa7osaNWRw/0hX6W/aPE18oG7t//zM3NVUxMDMHZbARnADCP3W7X+PHj9corr0iSbr75Zs2aNUtNmzY1raZz5845h08cP35cx44d04kTJ1yGFRQWFqq4uNhlmEJwcLAiIiJUr149l6VBgwZq2LCh6tWrV+VBEKgpqiqv8eVAAEC1VFJSomHDhmnu3LmSpCeeeELPPvus6T2uoaGhSklJUUpKiql1AKh6BGcAQLVjs9l01113acGCBQoMDNScOXP017/+1eyyAFzhCM4AgGrFMAyNGDFCCxYsUHBwsBYuXKg+ffqYXRYAcMttAED18swzz+jDDz9UYGCgPv30U0IzgGqD4AwAqDY++eQTPfPMM5KkGTNmqF+/fiZXBAD/h+AMAKgWdu/erREjRkiSHnnkEedjAKguCM4AANOdO3dOd9xxhwoKCtStWzdNnjzZ7JIAoByCMwDAdP/4xz+0Y8cOxcXF6eOPP1ZgIN9dB1D9EJwBAKb67rvvNG3aNEnSe++9p/j4eJMrAoCKEZwBAKYpKirSyJEjZRiG7rnnHt16661mlwQAF0RwBgCY5sUXX9Qvv/yi+Ph45221AaC6IjgDAEyxf/9+vfDCC5Kk1157TdHR0eYWBACXQHAGAJhi/PjxKioqUlpamv7yl7+YXQ4AXBLBGQBQ5b755ht99tlnslqtevXVV2WxWMwuCQAuieAMAKhShmFo/PjxkqR7771XqampJlcEAO4hOAMAqtSXX36p9evXKywsTE8//bTZ5QCA2wjOAIAqYxiGMyw/+OCDiouLM7cgAKgEgjMAoMp8+eWX+umnnxQWFqZHH33U7HIAoFIIzgCAKmEYhiZNmiTpfG9zbGysyRUBQOUQnAEAVWLp0qXasGGDwsPD6W0G4JcIzgAAnzMMQ88995wkepsB+C+CMwDA53788UetW7dOwcHBGjdunNnlAIBHCM4AAJ+bPn26JOmuu+5iJg0AfovgDADwqaNHj2rBggWSpP/6r/8yuRoA8BzBGQDgU2+99ZZKS0t18803q3Xr1maXAwAeIzgDAHzm7NmzevfddyVJDz/8sLnFAMBlIjgDAHxm7ty5ys7OVnJysv70pz+ZXQ4AXBaCMwDAJwzD0Ouvvy5JeuihhxQQEGByRQBweQjOAACf+Omnn7Rt2zaFhoZqxIgRZpcDAJeN4AwA8IlZs2ZJkv785z8rOjra3GIAwAsIzgAAryssLNTHH38sSbrnnntMrgYAvIPgDADwus8++0x5eXlq3LixunbtanY5AOAVBGcAgNc5hmkMHz5cViv/1ACoGfg0AwB41YEDB7R69WpZLBYNHTrU7HIAwGsIzgAAr5ozZ44k6Y9//KOuuuoqc4sBAC8iOAMAvKa0tNQZnPlSIICahuAMAPCatWvX6siRI4qKilK/fv3MLgcAvIrgDADwmnnz5kmSBg4cqNDQUJOrAQDvIjgDALyipKRE8+fPlyTdeeedJlcDAN5HcAYAeMWqVauUnZ2t+vXrq3v37maXAwBeR3AGAHiFY5jGHXfcocDAQJOrAQDvIzgDAC7buXPntHDhQkkM0wBQcxGcAQCXbdmyZcrLy1PDhg31H//xH2aXAwA+QXAGAFw2xzCNQYMGcYttADUWn24AgMtSUFCgf//735IYpgGgZiM4AwAuy7Jly1RYWKgmTZqobdu2ZpcDAD5DcAYAXBbHlwJvv/12WSwWk6sBAN8hOAMAPFZSUqLPP/9c0vngDAA1GcEZAOCxNWvWKDc3V7GxserQoYPZ5QCAT/ldcH7zzTeVnJys0NBQtW/fXuvXr79g2zlz5shisbgsoaGhVVgtANRsixYtkiT169dPAQEB5hYDAD7mV8H5k08+0dixYzVp0iRt2rRJrVq1Us+ePXXq1KkLbhMZGakTJ044l8OHD1dhxQBQc9ntdmdw7t+/v6m1AEBV8KvgPG3aNI0aNUrDhw/XtddeqxkzZig8PFyzZs264DYWi0Xx8fHOJS4urgorBoCaa8OGDTp+/Lhq166tHj16mF0OAPic3wTn4uJibdy4UWlpac51VqtVaWlpSk9Pv+B2+fn5atSokZKSktSvXz/9/PPPF32foqIi5eXluSwAgPIcs2n06dNHISEhJlcDAL7nN8E5KytLpaWl5XqM4+LilJmZWeE2zZs316xZs7R48WL9z//8j+x2uzp16qSjR49e8H0mT56sqKgo55KUlOTV4wCAmsIxTIPZNABcKfwmOHuiY8eOGjJkiFq3bq2uXbvqs88+U/369fXOO+9ccJsJEyYoNzfXuRw5cqQKKwYA/7Bnzx7t3r1bQUFB6tOnj9nlAECVCDS7AHfFxMQoICBAJ0+edFl/8uRJxcfHu7WPoKAgtWnTRvv27btgm5CQEP6XIwBcguMW2927d1dkZKTJ1QBA1fCbHufg4GC1bdtWq1atcq6z2+1atWqVOnbs6NY+SktLtX37diUkJPiqTAC4IixdulSSdOutt5pcCQBUHb/pcZaksWPHaujQobrxxhvVrl07vfbaayooKNDw4cMlSUOGDFGDBg00efJkSdKzzz6rDh06KCUlRadPn9bLL7+sw4cPa+TIkWYeBgD4tdzcXH333XeSxDANAFcUvwrOgwYN0q+//qqJEycqMzNTrVu31rJly5xfGMzIyJDV+n+d6L/99ptGjRqlzMxM1alTR23bttUPP/yga6+91qxDAAC/t2LFCtlsNjVv3lxNmzY1uxwAqDIWwzAMs4uozvLy8hQVFaXc3FzG8QGApHvuuUezZ8/W2LFj9corr5hdDgBUWV7zmzHOAADz2e125/hmhmkAuNIQnAEAbtu0aZNOnjypiIgIdenSxexyAKBKEZwBAG774osvJEl//OMfFRwcbHI1AFC1CM4AALcxTAPAlYzgDABwy6lTp/TTTz9Jknr37m1yNQBQ9QjOAAC3LFu2TIZhqE2bNkpMTDS7HACocgRnAIBbvvrqK0n0NgO4chGcAQCXZLfbncH5lltuMbkaADAHwRkAcElbt27Vr7/+qtq1a6tjx45mlwMApiA4AwAuydHb3L17d6ahA3DFIjgDAC5p+fLlkhimAeDKRnAGAFxUQUGB1q5dK0nq2bOnydUAgHkIzgCAi/rmm29UUlKi5ORkpaSkmF0OAJiG4AwAuKiywzQsFovJ1QCAeQjOAICLcnwxkGEaAK50BGcAwAVlZGRo9+7dslqt+sMf/mB2OQBgKoIzAOCCVqxYIUlq3769oqOjzS0GAExGcAYAXNDKlSslSX/84x9NrgQAzEdwBgBUyDAMff3115KkHj16mFwNAJiP4AwAqNDPP/+sU6dOKSwsTO3btze7HAAwHcEZAFAhR29zly5dFBISYnI1AGA+gjMAoEKO4MxsGgBwHsEZAFCOzWbTmjVrJBGcAcCB4AwAKGfz5s3Kzc1VVFSUbrjhBrPLAYBqgeAMACjHMUyjW7duCggIMLkaAKgeCM4AgHJWrVoliWnoAKAsgjMAwEVRUZHWrl0rifHNAFAWwRkA4OLHH3/U2bNnFRcXp2uvvdbscgCg2iA4AwBclJ2GzmKxmFwNAFQfBGcAgAvmbwaAihGcAQBO+fn5+vHHHyURnAHg9wjOAACntWvXymazKTk5WU2aNDG7HACoVgjOAAAnhmkAwIURnAEATo75mwnOAFAewRkAIEnKycnR5s2bJRGcAaAilxWci4uLtWfPHtlsNm/VAwAwyXfffSfDMNSiRQslJCSYXQ4AVDseBefCwkKNGDFC4eHhuu6665SRkSFJGj16tKZMmeLVAgEAVeObb76RJHXt2tXkSgCgevIoOE+YMEFbt27VmjVrFBoa6lyflpamTz75xGvFAQCqzpo1ayRJ3bp1M7UOAKiuAj3ZaNGiRfrkk0/UoUMHl7tKXXfdddq/f7/XigMAVI3Tp09ry5YtkuhxBoAL8ajH+ddff1VsbGy59QUFBdyeFQD80Nq1a2UYhpo1a8b4ZgC4AI+C84033qgvvvjC+dwRlt977z117NjRO5UBAKoM45sB4NI8GqrxwgsvqHfv3tq5c6dsNpumT5+unTt36ocffnB++AIA/AfjmwHg0jzqce7cubO2bNkim82mli1b6quvvlJsbKzS09PVtm1bb9cIAPChvLw8bdq0SRI9zgBwMR71OEtS06ZNNXPmTG/WAgAwwffffy+73a4mTZqoYcOGZpcDANWWRz3OmzZt0vbt253PFy9erP79++uJJ55QcXGx14oDAPge45sBwD0eBef77rtPv/zyiyTpwIEDGjRokMLDwzV//nyNHz/eqwUCAHzLEZwZ3wwAF+dRcP7ll1/UunVrSdL8+fPVtWtXzZ07V3PmzNH//u//erM+AIAP5efn66effpJEjzMAXIpHwdkwDNntdknSypUr1adPH0lSUlKSsrKyvFcdAMCnfvjhB5WWlqpRo0Zq1KiR2eUAQLXm8TzOzz//vP77v/9b33zzjfr27StJOnjwoOLi4rxaIADAdxjfDADu8yg4v/baa9q0aZMeeughPfnkk0pJSZEkLViwQJ06dfJqgQAA32F8MwC4z2IYhuGtnZ07d04BAQEKCgry1i5Nl5eXp6ioKOXm5ioyMtLscgDAawoLCxUdHa2SkhLt379fTZo0MbskAPBIVeU1j+dxlqTi4mKdOnXKOd7Z4aqrrrqsogAAvpeenq6SkhI1bNhQjRs3NrscAKj2PArOv/zyi0aMGKEffvjBZb1hGLJYLCotLfVKcQAA3yk7vtlisZhcDQBUfx4F5+HDhyswMFCff/65EhIS+MAFAD/EFwMBoHI8Cs5btmzRxo0b1aJFC2/XAwCoAufOndO6desk8cVAAHCXR7NqXHvttczXDAB+7Mcff1RRUZESEhKcMyMBAC7Oo+D84osvavz48VqzZo2ys7OVl5fnsgAAqjfGNwNA5Xk0VCMtLU2S1KNHD5f1fDkQAPwD45sBoPI8Cs6rV6/2dh0AgCpSVFSk9PR0SYxvBoDK8Cg4m9lD8eabb+rll19WZmamWrVqpTfeeEPt2rW7YPv58+frqaee0qFDh9SsWTO9+OKL6tOnTxVWDADVy/r163Xu3DnFxcWpefPmZpcDAH7DozHOknT69Gm98sorGjlypEaOHKlXX31Vubm53qytnE8++URjx47VpEmTtGnTJrVq1Uo9e/bUqVOnKmz/ww8/aPDgwRoxYoQ2b96s/v37q3///tqxY4dP6wSA6swxTOPmm29mfDMAVIJHt9zesGGDevbsqbCwMGdv708//aSzZ8/qq6++0g033OD1QiWpffv2uummm/Svf/1LkmS325WUlKTRo0fr8ccfL9d+0KBBKigo0Oeff+5c16FDB7Vu3VozZsxw6z0dt3A8fpxbbgOoGW677TatXv21Xnllmu677z6zywGAy5aXl6fERN/fctuj4NylSxelpKRo5syZCgw8P9rDZrNp5MiROnDggL799luvF1pcXKzw8HAtWLBA/fv3d64fOnSoTp8+rcWLF5fb5qqrrtLYsWP18MMPO9dNmjRJixYt0tatWyt8n6KiIhUVFTmf5+XlKSkpSVKuJIIzAABA9ZMnyffB2aOhGhs2bNBjjz3mDM2SFBgYqPHjx2vDhg1eK66srKwslZaWKi4uzmV9XFycMjMzK9wmMzOzUu0lafLkyYqKinIu50MzAAAArnQefTkwMjJSGRkZ5e4ceOTIEUVERHilMLNMmDBBY8eOdT539DgfPy4xUgOAv5s6daqefnqS/vSnfpo7d67Z5QCAV+TlSYmJvn8fj4LzoEGDNGLECE2dOlWdOnWSJH3//fd69NFHNXjwYK8W6BATE6OAgACdPHnSZf3JkycVHx9f4Tbx8fGVai9JISEhCgkJKbe+Vq3zCwD4s/T0lZIK1aNHBz7TANQYVXULEY+GakydOlUDBgzQkCFDlJycrOTkZA0bNkx//vOf9eKLL3q7RklScHCw2rZtq1WrVjnX2e12rVq1Sh07dqxwm44dO7q0l6QVK1ZcsD0A1GQ2m03ff/+9JOZvBgBPePTlQIfCwkLt379fktS0aVOFh4d7rbCKfPLJJxo6dKjeeecdtWvXTq+99po+/fRT7d69W3FxcRoyZIgaNGigyZMnSzo/HV3Xrl01ZcoU9e3bV/PmzdMLL7ygTZs2KTU11a33dMyq4evB5gDga+vXr1f79u1Vt25d/frrr7JaPZ6RFACqlarKax4N1XAIDw9XdHS087GvDRo0SL/++qsmTpyozMxMtW7dWsuWLXN+ATAjI8PlH4JOnTpp7ty5+sc//qEnnnhCzZo106JFi9wOzQBQkzjmb+7SpQuhGQA84FGPs81m0zPPPKPXX39d+fn5kqTatWtr9OjRmjRpkoKCgrxeqFnocQZQU/Tt21dLly7VtGnTNGbMGLPLAQCvqdY9zqNHj9Znn32ml156yTleOD09XU8//bSys7P19ttve7VIAMDlKS0t1dq1ayUxvhkAPOVRj3NUVJTmzZun3r17u6xfunSpBg8e7PNbb1clepwB1AQbN27UjTfeqKioKGVnZysgIMDskgDAa6oqr3k0yC0kJETJycnl1jdu3FjBwcGXWxMAwMvKjm8mNAOAZzwKzg899JCee+45l1tTFxUV6Z///KceeughrxUHAPCONWvWSGKYBgBcDo/GOG/evFmrVq1Sw4YN1apVK0nS1q1bVVxcrB49emjAgAHOtp999pl3KgUAeKS0tFTffvutJIIzAFwOj4JzdHS0Bg4c6LIuKSnJKwUBALxr69atznF/rVu3NrscAPBbHgXn2bNne7sOAICPOIZp3HzzzYxvBoDL4NEY57Nnz6qwsND5/PDhw3rttdf01Vdfea0wAIB3ML4ZALzDo+Dcr18/ffjhh5Kk06dPq127dnrllVfUr18/5nAGgGqE8c0A4D0eBedNmzapS5cukqQFCxYoPj5ehw8f1ocffqjXX3/dqwUCADzH+GYA8B6PgnNhYaEiIiIkSV999ZUGDBggq9WqDh066PDhw14tEADgudWrV0tifDMAeINHwTklJUWLFi3SkSNHtHz5ct1yyy2SpFOnTnF3PQCoRhjfDADe41Fwnjhxoh555BElJyerXbt26tixo6Tzvc9t2rTxaoEAAM8wvhkAvMtiGIbhyYaZmZk6ceKEWrVqJav1fP5ev369IiMj1aJFC68Waaaquvc5AHjbxo0bdeONNyoqKkrZ2dkM1QBQY1VVXvOox1mS4uPjFRERoRUrVujs2bOSpJtuuqlGhWYA8GfM3wwA3uVRcM7OzlaPHj109dVXq0+fPjpx4oQkacSIERo3bpxXCwQAeIbxzQDgXR4F5zFjxigoKEgZGRkKDw93rh80aJCWLVvmteIAAJ5hfDMAeJ9Ht9z+6quvtHz5cjVs2NBlfbNmzZiODgCqgS1btjjH/LVq1crscgCgRvCox7mgoMClp9khJydHISEhl10UAODyML4ZALzPo+DcpUsX5y23Jclischut+ull15S9+7dvVYcAMAzjG8GAO/zaKjGyy+/rD/84Q/asGGDiouLNX78eP3888/KycnR999/7+0aAQCVwPhmAPCNSgfnkpIS/f3vf9e///1vrVixQhEREcrPz9eAAQP04IMPKiEhwRd1AgDcxPhmAPCNSgfnoKAgbdu2TXXq1NGTTz7pi5oAAJeB8c0A4BsejXH+29/+pvfff9/btQAAvIDxzQDgGx6NcbbZbJo1a5ZWrlyptm3bqlatWi6vT5s2zSvFAQAqh/HNAOA7HgXnHTt26IYbbpAk/fLLLy6vWSyWy68KAOARxjcDgO94FJxXr17t7ToAAF7A+GYA8B2PxjgDAKqnr7/+WpLUtWtXkysBgJqH4AwANURJSYm++eYbSVJaWprJ1QBAzUNwBoAaYt26dSooKFBMTIxatmxpdjkAUOMQnAGghli5cqUkqUePHrJa+XgHAG/jkxUAaohVq1ZJYpgGAPgKwRkAaoAzZ87oxx9/lERwBgBfITgDQA3w7bffymazqUmTJkpOTja7HACokQjOAFADOMY309sMAL5DcAaAGoDxzQDgewRnAPBzmZmZ2r59uywWi7p37252OQBQYxGcAcDPOe4W2Lp1a8XExJhcDQDUXARnAPBzjG8GgKpBcAYAP2YYBsEZAKoIwRkA/Ni+fft05MgRBQcHq3PnzmaXAwA1GsEZAPyYo7e5U6dOCg8PN7kaAKjZCM4A4McYpgEAVYfgDAB+qrS0VKtXr5ZEcAaAqkBwBgA/tXnzZv3222+KiopS27ZtzS4HAGo8gjMA+Knly5dLkrp3767AwECTqwGAmo/gDAB+atmyZZKkXr16mVwJAFwZCM4A4IdOnz6t9PR0SVLPnj1NrgYArgwEZwDwQytXrlRpaalatGih5ORks8sBgCsCwRkA/JBjmEbv3r1NrgQArhwEZwDwM4ZhML4ZAExAcAYAP7Njxw4dO3ZMYWFhuvnmm80uBwCuGARnAPAzjt7m7t27KzQ01ORqAODKQXAGAD/z5ZdfSmKYBgBUNYIzAPiRM2fOaO3atZL4YiAAVDWCMwD4kdWrV6ukpERNmzZVSkqK2eUAwBWF4AwAfuTzzz+XxDANADADwRkA/ITdbncG5z/96U8mVwMAVx6CMwD4iY0bN+rEiROqXbu2unbtanY5AHDF8ZvgnJOTo7/+9a+KjIxUdHS0RowYofz8/Itu061bN1ksFpfl/vvvr6KKAcC7/v3vf0uSevbsqZCQEJOrAYArT6DZBbjrr3/9q06cOKEVK1aopKREw4cP17333qu5c+dedLtRo0bp2WefdT4PDw/3dakA4BNLliyRxDANADCLXwTnXbt2admyZfrpp5904403SpLeeOMN9enTR1OnTlViYuIFtw0PD1d8fHxVlQoAPpGRkaGtW7fKarWqT58+ZpcDAFckvxiqkZ6erujoaGdolqS0tDRZrVatW7fuott+9NFHiomJUWpqqiZMmKDCwsKLti8qKlJeXp7LAgBmcwzT6NSpk2JiYkyuBgCuTH7R45yZmanY2FiXdYGBgapbt64yMzMvuN1dd92lRo0aKTExUdu2bdNjjz2mPXv26LPPPrvgNpMnT9YzzzzjtdoBwBscwfm2224zuRIAuHKZGpwff/xxvfjiixdts2vXLo/3f++99zoft2zZUgkJCerRo4f279+vpk2bVrjNhAkTNHbsWOfzvLw8JSUleVwDAFyuM2fOaPXq1ZIIzgBgJlOD87hx4zRs2LCLtmnSpIni4+N16tQpl/U2m005OTmVGr/cvn17SdK+ffsuGJxDQkL4tjqAauWrr75ScXGxUlJS1KJFC7PLAYArlqnBuX79+qpfv/4l23Xs2FGnT5/Wxo0b1bZtW0nS119/Lbvd7gzD7tiyZYskKSEhwaN6AcAMixYtknS+t9lisZhbDABcwfziy4HXXHONevXqpVGjRmn9+vX6/vvv9dBDD+nOO+90zqhx7NgxtWjRQuvXr5ck7d+/X88995w2btyoQ4cOacmSJRoyZIhuvvlmXX/99WYeDgC4raioyDm+eeDAgSZXAwBXNr8IztL52TFatGihHj16qE+fPurcubPeffdd5+slJSXas2ePc9aM4OBgrVy5UrfccotatGihcePGaeDAgc5/gADAH6xatUq5ublKSEhQx44dzS4HAK5ofjGrhiTVrVv3ojc7SU5OlmEYzudJSUn65ptvqqI0APCZBQsWSDrf22y1+k1fBwDUSHwKA0A1VVJS4hzf/Oc//9ncYgAABGcAqK6+/vpr/fbbb4qNjVXnzp3NLgcArngEZwCophzDNAYMGKCAgACTqwEAEJwBoBqy2WxauHChJIZpAEB1QXAGgGrom2++UXZ2turVq6euXbuaXQ4AQARnAKiW5s+fL0m6/fbbFRjoNxMgAUCNRnAGgGqmuLjYGZz/8pe/mFwNAMCB4AwA1czSpUuVk5OjxMRE/eEPfzC7HADA/0dwBoBq5n/+538kSXfddRezaQBANUJwBoBq5LffftO///1vSdLdd99tcjUAgLIIzgBQjcyfP1/FxcVq2bKlrr/+erPLAQCUQXAGgGrEMUyD3mYAqH4IzgBQTRw6dEjfffedLBaL7rrrLrPLAQD8DsEZAKoJR29zjx491KBBA5OrAQD8HsEZAKoBu92uOXPmSGKYBgBUVwRnAKgGVq1apf379ysqKkoDBw40uxwAQAUIzgBQDcyYMUOSNGTIENWqVcvkagAAFSE4A4DJjh8/rsWLF0uS7rvvPpOrAQBcCMEZAEz23nvvqbS0VF26dNF1111ndjkAgAsgOAOAiWw2m2bOnClJuv/++02uBgBwMQRnADDR0qVLdfToUcXExPClQACo5gjOAGCit99+W5J0zz33KCQkxORqAAAXQ3AGAJPs2LFDy5Ytk8Vi0b333mt2OQCASyA4A4BJXnrpJUnSgAED1LRpU5OrAQBcCsEZAEyQkZGhjz/+WJL02GOPmVwNAMAdBGcAMMG0adNks9nUvXt33XTTTWaXAwBwA8EZAKpYdna2cwo6epsBwH8QnAGgir355psqLCxU69atdcstt5hdDgDATQRnAKhCZ86c0RtvvCFJGj9+vCwWi8kVAQDcRXAGgCo0ffp0ZWVlKSUlRXfccYfZ5QAAKoHgDABVJCcnR1OnTpUkPfvsswoMDDS5IgBAZRCcAaCKTJkyRbm5uWrZsqUGDRpkdjkAgEoiOANAFThw4ICmT58uSXrhhRdktfLxCwD+hk9uAKgC48ePV3FxsdLS0tS3b1+zywEAeIDgDAA+tmLFCv3v//6vrFarpk2bxkwaAOCnCM4A4EPnzp3Tf/7nf0qSHnzwQbVs2dLkigAAniI4A4APPf/889q3b58SExP1/PPPm10OAOAyEJwBwEd++uknTZkyRZL0+uuvKzIy0uSKAACXg+AMAD5QWFioIUOGqLS0VHfeeacGDhxodkkAgMtEcAYAH3jooYe0e/duJSQk6M033zS7HACAFxCcAcDLZs2apdmzZ8tqteqjjz5S3bp1zS4JAOAFBGcA8KLvvvtODzzwgKTzt9Xu3r27yRUBALyF4AwAXrJnzx7dfvvtKi4u1sCBAzVhwgSzSwIAeBHBGQC84MCBA+rRo4eys7N100036cMPP+S22gBQw/CpDgCXadeuXerWrZuOHTuma6+9Vl988YXCw8PNLgsA4GUEZwC4DGvWrFHnzp115MgRtWjRQitXrlT9+vXNLgsA4AOBZhcA4MpTVFSk3Nxc53LmzBnZ7XYZhiHDMCRJQUFBioqKUnR0tKKiohQZGanAwOrzkVVaWqpXXnlFTzzxhEpLS9W+fXt9/vnniomJMbs0AICPVJ9/hYAazmaz6ezZszp79qwKCwt17tw5Wa1WBQYGllvCwsIUGhoqi8VidtmVUlxcrOPHj+vIkSM6cuSIjh8/ruPHj+vEiRPOx8ePH1dhYaFH+4+Li1OjRo2UnJys5ORkNW7cWNddd51atmyp6Oho7x7MRezYsUP33XeffvjhB0nS3XffrRkzZjA8AwBqOIvh6N5BhfLy8hQVFaVRo0YpKChIkpw9YmVP3e/XVfSaJFmtVlmtVgUEBDh/ln1c9qdhGM5euAv9dKdNRdsEBASUC2tBQUEXfR4YGCi73a7S0lLZbDbn8vvnNptNxcXFKikpqXC52GslJSXOcxIYGOhSp+NxcHCwQkJCFBwc7Fwu9tydthaLxaW2sjUWFxfr3LlzzsDrCL+Xev7710pKSip17QUEBCgiIkIRERGKjIx0Po6IiHD2xDqWss/LPo6MjPTKF9RKSkqUk5OjrKwsZWVl6eTJk85wXHY5efKkKvOR4jiWiIgIBQQEyGKxOBdHr/Tp06d19uzZS+4rKSlJLVu21PXXX6+WLVuqVatWat68uVd7qXft2qWXXnpJH374oex2uyIiIjRt2jSNGDHC7/7IAYCaxJHXcnNzFRkZ6bP3IThfguMXAXiTo0fZMIxyf3R4k8ViUUREhEuorl27doV/qEly/oHgWPLy8pSVlaXTp0+7/Z7BwcFq2LChc0lMTHRZEhISFBMT4wzL7iguLtbp06d19OhRHTp0SIcPH9ahQ4e0b98+7dixQxkZGRVuFxISotTUVLVu3VqtWrVS69atde2116pu3bpuBV3DMHTw4EEtX75c8+bN07fffut87Y477tDUqVN11VVXuXdiAAA+Q3CuJhy/iMcee8zZMymp3M+K1v3+NUdvb2lpqUpLS529txd6bLVaZbFYvP7TYrHIbrfLZrOppKTEJbSVfV7Ra78fWvD7nmvH86CgIAUHBysoKKjccqH1jh5uxzkq25vt+Fm2F7ioqEjFxcXOpexzd19zPLbb7c66yv50PA4NDVVYWJjCwsIUHh7ufPz75+68drEhGI5jLyws1JkzZ3TmzBnl5eU5Hzue5+Xl6fTp087F0TNb9vG5c+e8+t+CxWJR3bp1FRMTo/r166thw4ZKSkoqt8TExFT5NGynT5/Wjh07tG3bNm3fvl3btm3Ttm3blJ+fX2H76OhopaSkqGHDhqpXr57q1q2rsLAw5/WQlZWlX375RXv27HH5o8Fisah///4aP368OnToUEVHBwC4FIJzNVFVvwjA28oOdSgbqAsKCmS3251/oDl+SnL5AyEsLEwRERGKiYlRTEyM6tSp43YPcXVgt9t14MABbdmyRVu3bnX+PHLkSKX2ExgYqE6dOqlv374aPHiwkpKSfFQxAMBTBOdqguAM1CyFhYU6cOCA9u/frxMnTig7O1s5OTk6d+6cc/x7ZGSkrr76al199dVq1qyZwsLCzC4bAHARVZXXmFUDwBUlPDxcqampSk1NNbsUAICf4QYoAAAAgBsIzgAAAIAbCM4AAACAGwjOAAAAgBsIzgAAAIAb/CY4//Of/1SnTp0UHh6u6Ohot7YxDEMTJ05UQkKCwsLClJaWpr179/q2UAAAANRIfhOci4uLdccdd+iBBx5we5uXXnpJr7/+umbMmKF169apVq1a6tmzp9fvqAYAAICaz+9ugDJnzhw9/PDDLrfBrYhhGEpMTNS4ceP0yCOPSJJyc3MVFxenOXPm6M4773Tr/bgBCgAAQPVWVXnNb3qcK+vgwYPKzMxUWlqac11UVJTat2+v9PT0C25XVFSkvLw8lwUAAACoscE5MzNTkhQXF+eyPi4uzvlaRSZPnqyoqCjnkpSU5NM6AQAA4B9MDc6PP/64LBbLRZfdu3dXaU0TJkxQbm6uczly5EiVvj8AAACqp0Az33zcuHEaNmzYRds0adLEo33Hx8dLkk6ePKmEhATn+pMnT6p169YX3C4kJEQhISEevScAAABqLlODc/369VW/fn2f7Ltx48aKj4/XqlWrnEE5Ly9P69atq9TMHAAAAIDkR2OcMzIytGXLFmVkZKi0tFRbtmzRli1blJ+f72zTokULLVy4UJJksVj08MMP6/nnn9eSJUu0fft2DRkyRImJierfv79JRwEAAAB/ZWqPc2VMnDhRH3zwgfN5mzZtJEmrV69Wt27dJEl79uxRbm6us8348eNVUFCge++9V6dPn1bnzp21bNkyhYaGVmntAAAA8H9+N49zVWMeZwAAgOqNeZwBAACAaoTgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4geAMAAAAuIHgDAAAALiB4AwAAAC4IdDsAqo7wzAkSXl5eSZXAgAAgIo4cpojt/kKwfkSsrOzJUlJSUkmVwIAAICLyc7OVlRUlM/2T3C+hLp160qSMjIyfPqLqEny8vKUlJSkI0eOKDIy0uxy/AbnrfI4Z57hvFUe58wznLfK45x5Jjc3V1dddZUzt/kKwfkSrNbzw8CjoqK4gCspMjKSc+YBzlvlcc48w3mrPM6ZZzhvlcc584wjt/ls/z7dOwAAAFBDEJwBAAAANxCcLyEkJESTJk1SSEiI2aX4Dc6ZZzhvlcc58wznrfI4Z57hvFUe58wzVXXeLIav5+0AAAAAagB6nAEAAAA3EJwBAAAANxCcAQAAADcQnAEAAAA3XPHB+Z///Kc6deqk8PBwRUdHu7WNYRiaOHGiEhISFBYWprS0NO3du9elTU5Ojv76178qMjJS0dHRGjFihPLz831wBOao7PEdOnRIFoulwmX+/PnOdhW9Pm/evKo4JJ/z5Jro1q1bufNx//33u7TJyMhQ3759FR4ertjYWD366KOy2Wy+PJQqVdnzlpOTo9GjR6t58+YKCwvTVVddpb///e/Kzc11aVeTrrU333xTycnJCg0NVfv27bV+/fqLtp8/f75atGih0NBQtWzZUkuXLnV53Z3PuJqgMudt5syZ6tKli+rUqaM6deooLS2tXPthw4aVu6Z69erl68OoUpU5Z3PmzCl3PkJDQ13acK2VV9HnvsViUd++fZ1tavq19u233+q2225TYmKiLBaLFi1adMlt1qxZoxtuuEEhISFKSUnRnDlzyrWp7GdlhYwr3MSJE41p06YZY8eONaKiotzaZsqUKUZUVJSxaNEiY+vWrcaf/vQno3HjxsbZs2edbXr16mW0atXK+PHHH43vvvvOSElJMQYPHuyjo6h6lT0+m81mnDhxwmV55plnjNq1axtnzpxxtpNkzJ4926Vd2fPqzzy5Jrp27WqMGjXK5Xzk5uY6X7fZbEZqaqqRlpZmbN682Vi6dKkRExNjTJgwwdeHU2Uqe962b99uDBgwwFiyZImxb98+Y9WqVUazZs2MgQMHurSrKdfavHnzjODgYGPWrFnGzz//bIwaNcqIjo42Tp48WWH777//3ggICDBeeuklY+fOncY//vEPIygoyNi+fbuzjTufcf6usuftrrvuMt58801j8+bNxq5du4xhw4YZUVFRxtGjR51thg4davTq1cvlmsrJyamqQ/K5yp6z2bNnG5GRkS7nIzMz06UN11p52dnZLudsx44dRkBAgDF79mxnm5p+rS1dutR48sknjc8++8yQZCxcuPCi7Q8cOGCEh4cbY8eONXbu3Gm88cYbRkBAgLFs2TJnm8r+Hi7kig/ODrNnz3YrONvtdiM+Pt54+eWXnetOnz5thISEGB9//LFhGIaxc+dOQ5Lx008/Odt8+eWXhsViMY4dO+b12quat46vdevWxj333OOyzp3/QPyRp+esa9euxn/9139d8PWlS5caVqvV5R+jt99+24iMjDSKioq8UruZvHWtffrpp0ZwcLBRUlLiXFdTrrV27doZDz74oPN5aWmpkZiYaEyePLnC9n/5y1+Mvn37uqxr3769cd999xmG4d5nXE1Q2fP2ezabzYiIiDA++OAD57qhQ4ca/fr183ap1UZlz9ml/l3lWnPvWnv11VeNiIgIIz8/37mupl9rZbnzWT1+/Hjjuuuuc1k3aNAgo2fPns7nl/t7cLjih2pU1sGDB5WZmam0tDTnuqioKLVv317p6emSpPT0dEVHR+vGG290tklLS5PVatW6deuqvGZv88bxbdy4UVu2bNGIESPKvfbggw8qJiZG7dq106xZs2TUgKnGL+ecffTRR4qJiVFqaqomTJigwsJCl/22bNlScXFxznU9e/ZUXl6efv75Z+8fSBXz1n9Lubm5ioyMVGBgoMt6f7/WiouLtXHjRpfPI6vVqrS0NOfn0e+lp6e7tJfOXzOO9u58xvk7T87b7xUWFqqkpER169Z1Wb9mzRrFxsaqefPmeuCBB5Sdne3V2s3i6TnLz89Xo0aNlJSUpH79+rl8LnGtuXeM77//vu68807VqlXLZX1NvdY8canPNW/8HhwCL90EZWVmZkqSS1BxPHe8lpmZqdjYWJfXAwMDVbduXWcbf+aN43v//fd1zTXXqFOnTi7rn332Wf3hD39QeHi4vvrqK/3nf/6n8vPz9fe//91r9ZvB03N21113qVGjRkpMTNS2bdv02GOPac+ePfrss8+c+63oWnS85u+8ca1lZWXpueee07333uuyviZca1lZWSotLa3wGti9e3eF21zomin7+eVYd6E2/s6T8/Z7jz32mBITE13+Ie7Vq5cGDBigxo0ba//+/XriiSfUu3dvpaenKyAgwKvHUNU8OWfNmzfXrFmzdP311ys3N1dTp05Vp06d9PPPP6thw4Zca25ca+vXr9eOHTv0/vvvu6yvydeaJy70uZaXl6ezZ8/qt99+u+z/5h1qZHB+/PHH9eKLL160za5du9SiRYsqqsg/uHveLtfZs2c1d+5cPfXUU+VeK7uuTZs2Kigo0Msvv1xtw4yvz1nZsNeyZUslJCSoR48e2r9/v5o2berxfs1WVddaXl6e+vbtq2uvvVZPP/20y2v+dq2h+pgyZYrmzZunNWvWuHzZ7c4773Q+btmypa6//no1bdpUa9asUY8ePcwo1VQdO3ZUx44dnc87deqka665Ru+8846ee+45EyvzH++//75atmypdu3auaznWjNPjQzO48aN07Bhwy7apkmTJh7tOz4+XpJ08uRJJSQkONefPHlSrVu3drY5deqUy3Y2m005OTnO7asjd8/b5R7fggULVFhYqCFDhlyybfv27fXcc8+pqKjI5/ef90RVnTOH9u3bS5L27dunpk2bKj4+vty3gk+ePClJV/y1dubMGfXq1UsRERFauHChgoKCLtq+ul9rFYmJiVFAQIDzd+5w8uTJC56f+Pj4i7Z35zPO33ly3hymTp2qKVOmaOXKlbr++usv2rZJkyaKiYnRvn37/D7MXM45cwgKClKbNm20b98+SVxrlzpvBQUFmjdvnp599tlLvk9NutY8caHPtcjISIWFhSkgIOCyr1+nSo2IrsEq++XAqVOnOtfl5uZW+OXADRs2ONssX768xn050NPj69q1a7kZDi7k+eefN+rUqeNxrdWFt66JtWvXGpKMrVu3Gobxf18OLPut4HfeeceIjIw0zp07570DMImn5y03N9fo0KGD0bVrV6OgoMCt9/LXa61du3bGQw895HxeWlpqNGjQ4KJfDrz11ltd1nXs2LHclwMv9hlXE1T2vBmGYbz44otGZGSkkZ6e7tZ7HDlyxLBYLMbixYsvu97qwJNzVpbNZjOaN29ujBkzxjAMrrVLnbfZs2cbISEhRlZW1iXfo6Zda2XJzS8HpqamuqwbPHhwuS8HXs7166ynUq1roMOHDxubN292To22efNmY/PmzS5TpDVv3tz47LPPnM+nTJliREdHG4sXLza2bdtm9OvXr8Lp6Nq0aWOsW7fOWLt2rdGsWbMaNx3dxY7v6NGjRvPmzY1169a5bLd3717DYrEYX375Zbl9LlmyxJg5c6axfft2Y+/evcZbb71lhIeHGxMnTvT58VSFyp6zffv2Gc8++6yxYcMG4+DBg8bixYuNJk2aGDfffLNzG8d0dLfccouxZcsWY9myZUb9+vVr3HR0lTlvubm5Rvv27Y2WLVsa+/btc5muyWazGYZRs661efPmGSEhIcacOXOMnTt3Gvfee68RHR3tnGnl7rvvNh5//HFn+++//94IDAw0pk6dauzatcuYNGlShdPRXeozzt9V9rxNmTLFCA4ONhYsWOByTTn+rThz5ozxyCOPGOnp6cbBgweNlStXGjfccIPRrFmzGvFHrGFU/pw988wzxvLly439+/cbGzduNO68804jNDTU+Pnnn51tuNbKnzeHzp07G4MGDSq3/kq41s6cOePMY5KMadOmGZs3bzYOHz5sGIZhPP7448bdd9/tbO+Yju7RRx81du3aZbz55psVTkd3sd+Du6744Dx06FBDUrll9erVzjb6//O9OtjtduOpp54y4uLijJCQEKNHjx7Gnj17XPabnZ1tDB482Khdu7YRGRlpDB8+3CWM+7tLHd/BgwfLnUfDMIwJEyYYSUlJRmlpabl9fvnll0br1q2N2rVrG7Vq1TJatWplzJgxo8K2/qiy5ywjI8O4+eabjbp16xohISFGSkqK8eijj7rM42wYhnHo0CGjd+/eRlhYmBETE2OMGzfOZdo1f1fZ87Z69eoK/5uWZBw8eNAwjJp3rb3xxhvGVVddZQQHBxvt2rUzfvzxR+drXbt2NYYOHerS/tNPPzWuvvpqIzg42LjuuuuML774wuV1dz7jaoLKnLdGjRpVeE1NmjTJMAzDKCwsNG655Rajfv36RlBQkNGoUSNj1KhRlf5HubqrzDl7+OGHnW3j4uKMPn36GJs2bXLZH9daxf+N7t6925BkfPXVV+X2dSVcaxf6HHecp6FDhxpdu3Ytt03r1q2N4OBgo0mTJi65zeFivwd3WQzDz+ZfAgAAAEzAPM4AAACAGwjOAAAAgBsIzgAAAIAbCM4AAACAGwjOAAAAgBsIzgAAAIAbCM4AAACAGwjOAAAAgBsIzgAAAIAbCM4AAACAGwjOAHCFWLZsmTp37qzo6GjVq1dPt956q/bv3292WQDgNwjOAHCFKCgo0NixY7VhwwatWrVKVqtVt99+u+x2u9mlAYBfsBiGYZhdBACg6mVlZal+/fravn27UlNTzS4HAKo9epwB4Aqxd+9eDR48WE2aNFFkZKSSk5MlSRkZGeYWBgB+ItDsAgAAVeO2225To0aNNHPmTCUmJsputys1NVXFxcVmlwYAfoHgDABXgOzsbO3Zs0czZ85Uly5dJElr1641uSoA8C8EZwC4AtSpU0f16tXTu+++q4SEBGVkZOjxxx83uywA8CuMcQaAK4DVatW8efO0ceNGpaamasyYMXr55ZfNLgsA/AqzagAAAABuoMcZAAAAcAPBGQAAAHADwRkAAABwA8EZAAAAcAPBGQAAAHADwRkAAABwA8EZAAAAcAPBGQAAAHADwRkAAABwA8EZAAAAcAPBGQAAAHADwRkAAABww/8DNWRJ9IiTAusAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from pyqsp.angle_sequence import Polynomial, QuantumSignalProcessingPhases\n",
    "\n",
    "poly = Polynomial(pcoefs)\n",
    "measurement = \"x\"\n",
    "ang_seq = QuantumSignalProcessingPhases(\n",
    "    poly, signal_operator=\"Wx\", method=\"laurent\", measurement=measurement\n",
    ")\n",
    "pyqsp.response.PlotQSPResponse(ang_seq, signal_operator=\"Wx\", measurement=measurement)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4835e5a8-5e24-44bf-a3ee-087675d6a60a",
   "metadata": {},
   "source": [
    "As $R(a)=-i*e^{i\\frac{\\pi}{4}Z}W(a)e^{i\\frac{\\pi}{4}Z}$ and we have odd number of rotations, we get an $i$ phase to our polynomial, so we get $Im(P(a))$ instead of the real part. So we will get the result in the $|1\\rangle$ state in the ancilla. However, we can fix it by adding $\\pi/2$ phase to the last or first rotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "de17bdc1-b753-40ca-b8f0-abe625eae462",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:36.305612Z",
     "iopub.status.busy": "2024-05-07T14:31:36.304368Z",
     "iopub.status.idle": "2024-05-07T14:31:36.310713Z",
     "shell.execute_reply": "2024-05-07T14:31:36.310015Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#### change the R(x) to W(x), as the phases are in the W(x) conventions\n",
    "phases = np.array(ang_seq)\n",
    "phases = phases - np.pi / 2\n",
    "phases[0] = phases[0] + np.pi / 4\n",
    "phases[-1] = phases[-1] + np.pi / 4 + np.pi / 2\n",
    "\n",
    "phases = (\n",
    "    -2 * phases\n",
    ")  # verify conventions. minus is due to exp(-i*phi*z) in qsvt in comparison to qsp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bee148a-9f0b-415f-a64a-173e2b23ceda",
   "metadata": {},
   "source": [
    "## Create the full QSVT Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "42952ee7-7fe8-4474-afe6-9f6de881fbd7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:36.315276Z",
     "iopub.status.busy": "2024-05-07T14:31:36.314111Z",
     "iopub.status.idle": "2024-05-07T14:31:36.321486Z",
     "shell.execute_reply": "2024-05-07T14:31:36.320762Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq import Constraints, create_model, synthesize\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def main(state: Output[QArray[QBit]], aux: Output[QBit]):\n",
    "    allocate(REG_SIZE, state)\n",
    "\n",
    "    hadamard_transform(state)\n",
    "    qsvt_fpaa(\n",
    "        phase_seq=list(phases),\n",
    "        arith_oracle=lambda state_reg, aux_reg: arith_equation(\n",
    "            state_reg[0 : REG_SIZES[\"a\"]],\n",
    "            state_reg[REG_SIZES[\"a\"] : REG_SIZES[\"a\"] + REG_SIZES[\"b\"]],\n",
    "            state_reg[REG_SIZES[\"a\"] + REG_SIZES[\"b\"] : REG_SIZE],\n",
    "            aux_reg,\n",
    "        ),\n",
    "        qbv_reg=state,\n",
    "        qbv_aux=aux,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6493595c-d92d-4470-865c-367b1ba0e662",
   "metadata": {},
   "source": [
    "## Synthesis and execution on a simulator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1de3b81-2152-4ecf-9cd4-97be748561f9",
   "metadata": {},
   "source": [
    "We will use Classiq's synthesis engine to translate the model to a quantum circuit, and execute on Classiq's simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8591b0f2-d81a-4fbd-8f95-e6f0017f7b6e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:36.326397Z",
     "iopub.status.busy": "2024-05-07T14:31:36.324931Z",
     "iopub.status.idle": "2024-05-07T14:31:36.516252Z",
     "shell.execute_reply": "2024-05-07T14:31:36.515336Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import (\n",
    "    create_model,\n",
    "    execute,\n",
    "    set_execution_preferences,\n",
    "    show,\n",
    "    synthesize,\n",
    "    write_qmod,\n",
    ")\n",
    "from classiq.execution import (\n",
    "    ClassiqBackendPreferences,\n",
    "    ClassiqSimulatorBackendNames,\n",
    "    ExecutionPreferences,\n",
    ")\n",
    "\n",
    "# convert the functions to a qmod model\n",
    "constraints = Constraints(max_width=12)\n",
    "qmod = create_model(main, constraints)\n",
    "\n",
    "NUM_SHOTS = 1000\n",
    "# we will want to execute this qmod on a state-vector simulator:\n",
    "execution_preferences = ExecutionPreferences(\n",
    "    num_shots=NUM_SHOTS,\n",
    "    backend_preferences=ClassiqBackendPreferences(\n",
    "        backend_name=ClassiqSimulatorBackendNames.SIMULATOR\n",
    "    ),\n",
    ")\n",
    "qmod = set_execution_preferences(qmod, execution_preferences)\n",
    "\n",
    "write_qmod(qmod, \"qsvt_fixed_point_amplitude_amplification\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "69c827b2-0330-40d3-af08-1d669e8c43a5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:31:36.519850Z",
     "iopub.status.busy": "2024-05-07T14:31:36.519627Z",
     "iopub.status.idle": "2024-05-07T14:32:52.629895Z",
     "shell.execute_reply": "2024-05-07T14:32:52.629168Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening: https://platform.classiq.io/circuit/52d4df69-abbb-4410-847e-ea2e8fcb3cf4?version=0.41.0.dev39%2B79c8fd0855\n"
     ]
    }
   ],
   "source": [
    "qprog = synthesize(qmod)\n",
    "show(qprog)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "762aec30-4746-4632-b3b4-c104cab731a2",
   "metadata": {},
   "source": [
    "Execute the circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4f0678db-c669-420a-8e17-b071e4b35173",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:32:52.634302Z",
     "iopub.status.busy": "2024-05-07T14:32:52.632934Z",
     "iopub.status.idle": "2024-05-07T14:33:07.266667Z",
     "shell.execute_reply": "2024-05-07T14:33:07.265833Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq.execution import ExecutionDetails\n",
    "\n",
    "raw_results = execute(qprog).result()\n",
    "results = raw_results[0].value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "93e7a0d5-ca40-4ffa-936c-da8a750408c2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:33:07.271437Z",
     "iopub.status.busy": "2024-05-07T14:33:07.270773Z",
     "iopub.status.idle": "2024-05-07T14:33:07.282566Z",
     "shell.execute_reply": "2024-05-07T14:33:07.281885Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a: 3, b: 0, c: 5, aux: 0, equation_result: True, counts=435\n",
      "a: 2, b: 1, c: 4, aux: 0, equation_result: True, counts=393\n",
      "Measured good shots: 828\n"
     ]
    }
   ],
   "source": [
    "def to_int(x):\n",
    "    return int(x[::-1], 2)\n",
    "\n",
    "\n",
    "def equation(a, b, c):\n",
    "    return eval(EQUATION)\n",
    "\n",
    "\n",
    "measured_good_shots = 0\n",
    "\n",
    "for keys, counts in results.counts_of_multiple_outputs([\"state\", \"aux\"]).items():\n",
    "    a, b, c, aux = (\n",
    "        to_int(keys[0][: REG_SIZES[\"a\"]]),\n",
    "        to_int(keys[0][REG_SIZES[\"a\"] : REG_SIZES[\"a\"] + REG_SIZES[\"b\"]]),\n",
    "        to_int(keys[0][-REG_SIZES[\"c\"] :]),\n",
    "        to_int(keys[1]),\n",
    "    )\n",
    "    if equation(a, b, c) and (aux == 0):\n",
    "        print(\n",
    "            f\"a: {a}, b: {b}, c: {c}, aux: {aux}, equation_result: {equation(a, b, c)}, counts={counts}\"\n",
    "        )\n",
    "        measured_good_shots += counts\n",
    "\n",
    "print(\"Measured good shots:\", measured_good_shots)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddb24bd0-ae8d-4b2e-9d4b-fc6aec0cd298",
   "metadata": {},
   "source": [
    "What do we expect?\n",
    "\n",
    "We need to subtitue the amplitude of $|s\\rangle\\langle t|$ in $P(x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d4387720-27dc-45ac-9938-98c9fceda58f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:33:07.287126Z",
     "iopub.status.busy": "2024-05-07T14:33:07.285952Z",
     "iopub.status.idle": "2024-05-07T14:33:07.292662Z",
     "shell.execute_reply": "2024-05-07T14:33:07.292006Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected good shots: 832.9873728406709\n"
     ]
    }
   ],
   "source": [
    "p_good_shot = poly(np.sqrt(2 / 2**6)) ** 2\n",
    "print(\"Expected good shots:\", NUM_SHOTS * p_good_shot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0af1b8bf-2e53-4c4b-bb31-27b7463bc245",
   "metadata": {},
   "source": [
    "Indeed, we got the expected result according to the polynomial we created with the QSVT sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9d366c6f-79fa-4675-9b8e-c852fff436f6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:33:07.297231Z",
     "iopub.status.busy": "2024-05-07T14:33:07.296014Z",
     "iopub.status.idle": "2024-05-07T14:33:07.490186Z",
     "shell.execute_reply": "2024-05-07T14:33:07.489502Z"
    }
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "\n",
    "assert np.isclose(\n",
    "    measured_good_shots,\n",
    "    NUM_SHOTS * p_good_shot,\n",
    "    atol=5 * scipy.stats.binom.std(NUM_SHOTS, p_good_shot),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "982cd6a4-7abc-4f87-a00d-12e5445cc94b",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<a id='grand'>[1]</a>: [Martyn JM, Rossi ZM, Tan AK, Chuang IL. Grand unification of quantum algorithms. PRX Quantum. 2021 Dec 3;2(4):040203.](https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.040203)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
